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CHAPTER I 


INTRODUCTION AND STATEMENT OF PROBLEM 

Atmospheric turbulence and mean winds have influ- 
enced the design of the control, stability and guidance of 
aircraft for many years. Recently, the study of atmospheric 
flow around buildings and other man-made obstructions has 
attained major importance in the design and operational 
procedures of helicopter and V/STOL aircraft operating in 
large metropolitan areas. Such low speed aircraft flying 
near buildings may experience fields of induced vorticity, 
zones of recirculation and regions of large fluctuations 
which can make takeoff and landing extremely hazardous. 

Civil engineers are also concerned with unsteady 
atmospheric flow which induces dynamic loads on structures. 
Pollution control and monitoring are strongly influenced by 
the motion of the atmosphere over buildings. Many other 
areas of current technology also require an understanding of 
atmospheric flows in the vicinity of buildings. Since full- 
scale experimental studies of building flow phenomena are 
costly and exact simulation in wind tunnels remains to be 
verified the need for analytical methods to predict local 
atmospheric motion influenced by surface obstructions is 
thus obvious. The purpose of this study is to develop a 
meaningful mathematical model of flow over a single two- 
dimensional block geometry. The problem concerned herein is 



the wind field over a rectangular block on a flat plane 
surface. Undisturbed neutrally stable atmospheric wind 
perpendicular to the axis of the building is assumed far up- 
stream and far downstream of the obstacle, see Figure 1. 

To provide a better understanding of fluid mechanics 
of the problem, the flow regions associated with bluff body 
flow are qualitatively discussed in the following. Frost [1] 
has compiled an extensive survey of flow properties around 
man-made surface obstructions to the wind. The distorted 
shear flows approaching and pa's sing over a block geometry 
can be divided into three basic flow regions. Figure 2; 

(1) the displacement zone, (2) the upstream separation 
bubble or downwash zone, and (3) the wake zone which includes 
the rear separation bubble or cavity zone. The effect of 
shear in the approaching flow generates a vortex near the 
ground upstream of the block, creating a downwash on the 
front face. This can be interpreted as the piling up of 
vortex lines swept in by the incident flow. Also, shear 
generates a swirling flow in the wake or cavity zone, which 
can be explained as the accumulation of vortex lines in the 
wake zone by the passing flow. When the flow passes over 
the abrupt change in geometry at the corner of the block 
surface, separation from the sharp edge occurs. On the for- 
ward surface of the block the strong pressure gradient above 
the region of the separated eddy accelerates the fluid in 

^Numbers in brackets refer to similarly niimbered 
references in the Bibliography. 
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Figure 1. Description of flow region considered. 
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Definition of flow zones near a sharp-edged 

[ 1 ] . 



the thin surface boundary layer toward the sharp leading 
edge where, due to its vertically directed momentum, it is 
forced to separate. The separated flow forms a shear layer 
of low static pressure and high vorticity which is bent down- 
stream through interaction with the transverse main flow 
forming an essential shell which reattaches some distance 
downstream. The intensity of turbulence in the incident flow 
tends to force the wake to reattach closer to the back-side 
of the block and to thicken the shear layer which bounds the 
wake zone. An adverse pressure gradient and viscosity are 
two necessary conditions for flow separation to occur on the 
upstream face [2] . Downstream of separation from the corner 
momentum in the separated layer diffuses into the wake and 
into the quasi-potential flow outside the wake setting the 
wake fluid into motion and smoothing out the sharp velocity 
discontinuity. Beyond the point of reattachment this diffu- 
sion gradually thickens the shear layer until the inner flow 
is blended with the outer flow forming a new and thicker 
boundary layer. At the forward separation point and at the 
reattachment point of the wake where the normal velocity 
gradient at the surface is zero, the shear stress in laminar 
flow is reduced to zero. The longitudinal position on the 
surface where the shear stress becomes zero is identified as 
the separation point or reattachment point. Experimentally, 
the above description of zero stress which correlates with 
the onset of reverse flow is only true for a steady two- 
dimensional, laminar incompressible flow. For turbulent flow 
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the zero stress position fluctuates widely. However, 
ensemble time-average quantities are generally employed such 
that the mean normal velocity gradient is zero, i.e. , the 
separation point and the reattachment point are assimed to 
coincide with the condition of zero normal gradient in the 
mean velocity. 

Wakes generated by obstacles are characterized by 
increased turbulence, a mean velocity defect, and in certain 
situations by organized, discrete standing vortices. One of 
the first efforts to determine the characteristics of full- 
scale building wake was made by Colmer [3] . Colmer drew the 
conclusion for his study that the mean velocity deficit had 
completely decayed by 14 obstacle heights downstream and 
turbulence intensity excess, although small, was still 
apparent at 14 heights downstream. Frost and Shahabi [5] 
conducted a field test of wind passing over a simulated 
rectangular block building and observed that the average 
midpoint of reattachment along the centerline of the wake 
was 12.5 ± 2.34 building heights. Typical values measured 
for two-dimensional bluff bodies in wind tunnels range from 
13 to 17 building heights [5] . Lemberg [6] noticed that the 
mean velocity wake had effectively decayed in 12 building 
heights and the turbulence wake extended to 50 heights behind 
a cubical obstacle on a plane surface in the wind tunnel. In 
general, three-dimensional wakes are smaller than two- 
dimensional wakes. It is reasonable to expect an obstacle 
having a height which is a significant fraction of the 
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boundary layer will create a larger relative disturbance to 
the boundary layer than the disturbance due to a small 
obstacle. Observations reported in [4] show that the mean 
velocity decay rate in the wake is independent of model scale 
size, but the magnitude of the disturbance is, in a non- 
dimensional sense, greater for the larger model and hence its 
wake persists farther downwind (Figure 3). Hunt [8] devel- 
oped a theory for the wake behind a cubic block geometry and 
compared it with Counihan's [8] experimental results. 

Initial comparison gave good results but insufficient experi- 
mental data was available to test the breadth of the theory. 

At present, no well-developed simple theories are 
available to describe the wakes behind two-dimensional obsta- 
cles. Townsend [9] investigated the effects of a change in 
wall boundary conditions on an equilibrium layer. In 
practice his analysis is only valid sufficiently far down- 
stream of the obstacles that no new length scale is required. 

Bradshaw and Wong [12] predicted that the length 
scale of turbulence increased rapidly above the local 
equilibrium value with increasing distance normal to the 
surface which evidently arises from the bifurcation of the 
mixing layer at reattachment. Thus Townsend's theory is 
valid for an ordinary boundary layer but not for a reattach- 
ing shear layer. Upstream from reattachment, however, the 
region of strong shear near the mean separation streamline 
is known to be reasonably described by the theory of a 
turbulent mixing layer [10] . However, Tani, et al. [27] 
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O Model 1 (z/H = 1.08) (H/6 = 0.05) 

A Model 2 (z/H = 0.94) (H/6 = 0.107) 

a Model 3 (z/H = 0.8) (H/6 = 0.25) 

• Colmer (z/H = 0.9) (H/6 = 0.017) 


Figure 3. Decay of mean velocity deficit along the wake 
center line for values of H/6 [4] . 




and Mueller and Robertson [7] measured shear stress in the 
free shear layer much higher than that found in a plane 
mixing layer. The reason for higher stress is that the 
effective velocity difference across the shear layer is more 
than the free stream velocity, due to the reversed flow 

in the separated region. Existing computational methods 
based on turbulent boundary layer concepts are not applicable 
to problems in which recirculation and separation occur. 

Since the topic of this study is directly concerned with the 
flow phenomena around bluff body surface obstructions to the 
wind and are complicated by separation and recirculation, 
the complete equations of motion must be considered if a 
realistic flow field is to be computed. Practical methods, 
however, for solving such equations are limited to numerical 
approaches which so far have assumed nonturbulent or turbu- 
lent flow based on a constant effective viscosity. More 
realistic turbulence models for closing the system of 
governing equations when Reynolds averaging has been carried 
out are in the development stages [21]. Considering 
ensemble averaged momentum equations for turbulent incom- 
pressible flow, one can have 


_ 3u^ 


u 


3 3x- 


_1 3P_ ^ 
P 


V 


3 ^Uj^ 


3x . 3x . 
1 3 


3 

3xj 



where over-bars denote the ensemble average and superscripts 
” denote fluctuation properties. The above equation can 
be readily integrated if certain assumptions on the Reynolds 
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stress, puTuT, can be made to form a closed set of equa- 
tions. Usually, Reynolds shear stress is expressed by the 
product of an eddy viscosity and a velocity gradient as 
follows : 


pu|u^ = 


3x . 
3 


3u. 


3x . 
1 


The eddy viscosity is a feature of turbulence so that it 
is a property of the flow field rather than that of the 
state of the fluid as is mole,cular viscosity. The value of 
eddy viscosity will thus vary from point to point in the 
flow, being largely determined by the structure and history 
of turbulence. Models for expressing the turbulent or eddy 
viscosity in terms of known or calculable quantities of the 
flow field are well described by Launder and Spalding [21] . 
Prandtl's mixing length model allows realistic prediction 
in boundary layer flows, but takes no account of the process 
of convection or diffusion of turbulence. The mixing length 
model takes 


y 


t 




3z 


so this means that the turbulent viscosity vanishes whenever 
the mean velocity gradient is zero. For the problems 
related to a flow downstream from a step, the mean velocity 
gradient is zero at the reattachment point; so if the mixing 
length hypothesis is adopted, the value of y^ is also zero 
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there. However, the level of velocity fluctuations in the 
neighborhood of reattachment point has been shown to be 
high [21] . Generally, the local level of velocity fluctua- 
tions is determined, not only by events at the point in 
question but also by influences which have originated some 
distance upstream. In a recirculating flow, downstream 
happenings can also influence the value of the turbulence 
at any point [21] . It is important to take account of the 
transport characters of turbulence in the problems for which 
the mixing length model is not well equipped. To obtain 
the realistic predictions of the recirculating flow situa- 
tion, the convective transport characters of turbulence 
have to be considered. A comparison of the performance of 
the different models on free shear flows [13, 20] indicates 
that the calculation of a length scale, £, from a transport 
equation gives more correct predictions over a wider range 
of flows than is achieved with models embodying an algebraic 
expression for £. Prandtl and Kolmogorov develop indepen- 
dently a turbulent model which formulates 
where both the turbulence kinetic energy, K, and the 
turbulence length scale, £, are calculated from transport 
equations, referred to as a two-equation model. This 
model is used throughout this study. Gosman, et al. [11] 
have written a complete volume on heat and mass transfer 
in two-dimensional recirculating flows. Their method has 
been applied to turbulent flows over a bluff surface 
obstruction [13] and over a backward facing step [14] . The 
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program developed in [13] is adapted in this study to the 
investigation of the flow over a bluff block. 

The numerical study reported in this investigation 
is being simultaneously supported by a field study [5] and 
the wind tunnel studies [4, 28,29] . The field study [5] is 
being conducted in the; NASA Marshall Space Flight Center, 
Space Sciences Laboratory, Atmospheric Sciences Division 
Atmospheric Boundary Layer Facility (ABLF) , located in 
Huntsville, Alabama. The tower array, located in a large 
open field, is schematically shown in Figure 4. There is a 
line of trees approximately 122 m upstream from the simulated 
block which are not shown. The existence of a wake from 
the tree line is clearly evident in the results of the field 
study [5] . The results from Reference [5] are summarized 
as follows: 

1. Analysis of the mean longitudinal velocity at 
building level shows an overshoot at tower 
number 3 . 

2. The smoke study shows that the mean extent of 
the wake is 12.5 +2.3 building heights under 
the conditions of approximately a 16 mph wind 
speed at the 20m level. 

3. The smoke study also indicates that the wake 
extends upward to a height of approximately 1.5h 
and 2. Oh and the separation at the upstream 
face of the building extends forward about 0.9h. 

4. The rms values of the velocity components at the 
3 m level were influenced by the building but at 
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Note : 


Measurements in meters 88 

Tower 3^ instrument level 4 was 
moved to the 9 meter level 



Figure 4. Tower arrangements carried out in field studies [5] 



the 12 m level this influence was not apparent, 
indicating that the disturbulence from the 
building did not extend to that height. 

5. The comparison of turbulence energy spectra and 
the ratio of velocity deviations, a^/a^/a^, to 
those of neutrally stable atmosphere over 
homogeneous uniform terrain reported in the 
literature show the influence of the tree wake 
to be evident. 

The wind tunnel tests carried out in Reference [4] 
indicate the tree line in the front of the simulated block 
gives the flow approaching the building a higher turbulence 
level. Figure 5, and a higher exponent in the power-law 
velocity profile (from 0.25 without the trees to 0.35 with 
trees). Figure 6. 


The power law velocity profile is described by 


u 


Uref 


'im 


^ref J 


is taken as 


which is the same as the logarithmic law if m 
1 


iln(Zj.ef/zo) * 

The data from Reference [4] which did not model the 
trees in most of the measurements made is directly comparable 
to the results of this numerical study if three-dimensional 
effects can be neglected. These effects will be discussed 
and comparisons will be made in Chapter IV. 

The governing equations to be solved are the Navier- 
Stokes equations with the Prandtl-Kolmogorov hypothesis for 
turbulent viscosity which requires simultaneous solution of 
two additional transport equations, the turbulence kinetic 
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o|n 





Figure 6. Vertical profile of longitudinal mean velocity in 
the approach flow [4] . 



energy and the turbulence length scale equations. The 
method of solution chosen is to eliminate the pressure terms 
involved in the equations of motion by employing the vor- 
ticity and the stream function as the dependent variable 
[11] . The solution is therefore obtained without any 
information being required about the pressure gradients. 
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CHAPTER II 


MATHEMATIC MODEL FOR ANALYSIS OF ATMOSPHERIC 
FLOW OVER A TWO-DIMENSIONAL 
RECTANGULAR BLOCK 

The complete two-dimensional Navier-Stokes equations 
of motion are applied to an atmospheric flow over a surface 
obstruction simulated by a sharp-edged rectangular block 
geometry as shown in Figure 1, page 3. Being a turbulent 
flow field, the ensemble time-average equations are used. 

The atmosphere is assumed to be neutrally stable and 
Coriolis effects induced by the rotation of the earth are 
considered negligible in the lower regions of the atmosphere 
to which the present problem is confined. 

1. MEAN VELOCITY PROFILE OF APPROACHING WIND 


For a neutral atmosphere, experimental evidence 
[15, 16] confirms the mean wind velocity in the region near 
the ground as being described by a logarithmic wind profile. 
The logarithmic wind profile is thus [17] 


u 



z + 2 , 


( 1 ) 


where Zq represents the surface roughness and k is von 
Karman's constant. Observed wind profiles up to 150 m over 
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reasonably level and uniformly rough terrain under neutrally 
stable conditions obey this law very well [18] . 

2. THE TURBULENT FLUXES NEAR THE GROUND 

Associated with the mean velocity gradient is a 
shear stress arising both from molecular or viscous momentum 
transport and from turbulence momentum transport. The 
turbulence momentum transport is generally much larger than 
the molecular transport. In the unstable atmospheric layer, 
heat transfer from the earth to the ambient air or thermal 
convection enhance intensity of the turbulence in the 
boundary. Conversely, heat transfer to the earth from the 
air decreases or even suppresses the turbulence in the 
boundary. In the neutral boundary layer, however, the 
effect of thermal convection is negligible and the mechanism 
of turbulence is solely due to the effects of shear. At a 
height z, small compared with the boundary layer thickness, 
the stress varies little from its surface value which is 
traditionally denoted by pu*. Reference [19] shows that: 


-pu'w' 


1 


f) 


pu| ~ pu^, 


z << 6 


(2) 


Equation 2 suggests that at height < 0.16, the shear stress 
is within 10% of its surface value. This has led to the 
concept of a "constant shear stress" layer near the earth's 
surface. Outdoor measurements have provided data which 
supports the concept of a constant-flux layer [19], Figure 7. 
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This stress picture is surely an idealized model since it is 
restricted to steady-state and horizontally homogeneous 
velocity conditions. However, for high wind speeds which 
are of primary importance to this study the concept of a 
neutrally, stable horizontally homogeneous and statistically 
steady boundary layer is meaningful. 

3. GOVERNING EQUATIONS 


The governing ensemble averaged turbulent flow 
equations for statistically steady incompressible two- 
dimensional flow can be written as follows: 

Continuity Equation: 


8u , ^ 

3x 3z 


0 


(3) 


Momentum Equations ; 


(1) X direction: 
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(4) 


(5) 


where the overbars denote ensemble averaged quantities. The 
"effective viscosity," is defined by 

^eff = ^ + ’'t 
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which is composed of molecular viscosity, y , and of turbu- 
lent viscosity, 

For two-dimensional flow, the stream function, 
and vorticity, u, can be introduced by the relationships 


1 

p 9z 




(7) 


p 9x 


( 8 ) 


— 9 w 9 u 

^ ~ 9x 9z 


(9) 


An equation known as the stream function equation is derived 
from Equations 7 through 9 

V^ip = -po) (10) 

Introducing Equations 7 through 9 into the momentum equa- 
tions and then differentiating the x-momentum equation with 
respect to z, the z-momentum equation with respect to x, and 
subtracting one from the other, an equation which is 
generally known as the vorticity transport equation is 
obtained. After some arranging, the governing equation for 
u) becomes 


9 

f- 

_9_ 


9x 


9z 






0 

( 11 ) 


where 
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( 12 ) 


Equation 11 reduces the two momentum equations in three 
variables u, w, p, to one equation with two variables ip, u. 
A turbulence closure problem arises, however, because 
Equation 11 contains three unknowns, ip, w, and and 

only two equations (10 and 11) are available for solution. 

A turbulence model for is therefore needed and its 

development is discussed in Section IV. 

Both of the new governing equations (10 and 11) can 
be expressed in the following general elliptical equation 
form (Equation 11) . 


3x 




dip 

3z 


_ 3 ^ 

3z 


^ 3x 


3 

3x 


b ^(ccO) 


3 

3z 


b j^(c<p) 


+ d = 0 
(13) 


where the coefficient functions, a, b, c, and d, correspond- 
ing to the dependent variable <() are given in Table I. 


4 . TURBULENCE MODEL 


The turbulent viscosity concept provides a framework 
for constructing a turbulent model, but there remains the 
task of formulating in terms of known or calculable 
quantities of the mean flow. 

In recirculating flows the convective transport of 5- 
is expected to be large. The investigations [13, 14, 21] 
indicated that the two-equation model adequately accounts 
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TABLE 1 



COEFFICIENT 

FUNCTIONS 

OF EQUATION 13 


(f> 

a 

b 

c 

d 


0 

1/P 

1 

-0) 

(A) 

1 

1 

^eff 

-^0, 

K 

1 

^K,eff 

1 

“^K 

£ 

1 

^a,eff 

1 
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for the convective character of turbulence in recirculating 
flows. The principal idea of the two-equation model is that 
the main local feature of a turbulent fluid can be repre- 
sented by just two quantities — the kinetic energy of the 
fluctuating motion, K, and the length scale of the turbu- 
lence, a. These two quantities are determined from 
appropriately determined transport equations. 

Kolmogorov and Prandtl proposed a differential 
equation for K assuming that is dependent on the level of 

turbulence of the fluid and the length scale as shown in 
equation as reproduced here for convenience. 

y . = C • £ (14) 

■^t yo® 

is an asymptotic coefficient to be discussed later. The 
level of turbulence is defined as the mean kinetic energy of 
the velocity fluctuation, or simply the kinetic energy of 
turbulence . 

K = i(u'2 + v'2 + w'2) (15) 

Combining Equations 6 and 14 and introducing a coefficient 
function 



one obtains 



(16) 


( 17 ) 
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An outline of the derivation of the transport equation for K 
is given by Wolfshtein [22] 


- 9K 

pu . 


j 9x . 9x . [ K, ef f 9x . 


9K 


+ S 


K 


(18) 


where is an exchange coefficient for the turbulence 

energy defined in terms of a Schmidt number based on the 
effective viscosity. 


^K,eff 


’^eff 
'^K, ef f 


(19) 


which represents "sources" and "sinks" of turbulence, is 
denoted as the following. 




8u . 
1 

9x . 
3 




9u . 


9x . 
1 


9u . 

^ 


K 


3 / 2 


( 20 ) 


is a coefficient function to be determined empirically. 

Rotta (1951) derived a differential equation for 
the length scale from the time-dependent motion. Gosman, 
et al. [11] reformulated the equation into the common 
elliptical form of their transport equation. The equation 
for I is thus: 


- 9£ 


9 

9x . 
3 


, 9£ 

5. , ef f 9x . J 
3^ 


+ S, 


( 21 ) 


where F^^ is an exchange coefficient for the length scale 
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in terms of a Sclimidt number based on the effective vis- 
cosity 


^il,eff " 


^ef f 


( 22 ) 


is a combined source and sink term for the turbulent 
length scale; 


= "^t 


■ah, 

ax. 

j 


auj au, 


ax . ax . 


I • s + p'' 


1/2 


(23) 


Again C_ and C_ are coefficient functions to be determined. 

x3 o 

Here, both the K and i transport equations are taken 
as obeying the general form of the elliptical equation, 
Equation 13, and the appropriate coefficient are also given 
in Table I, page 24. 


5. COEFFICIENT FUNCTIONS 


It has been assumed that the properties of local 
turbulence can be characterized by the kinetic energy of the 
fluctuating motion, K, and the length-scale of turbulence, 

Z. From these transport properties and from the density and 
viscosity of the fluid, a "Reynolds, number of turbulence," 
R^, is defined as 



Together with Equation 16, one can obtain 
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It is expected that when is very small, turbulence is 
negligible and the function tends to 1/R^; when R^ is 
very large, molecular viscosity has no significant effect 
and the function tends asymptotically to a constant 
value. An asymptotic value of for recirculating flows is 
determined in this report. 

The function is a measure of the dissipation rate 
per unit volume of turbulence kinetic energy. The coeffi- 
cient function is argued to be a function of R^ having 
similar behavior as that predicted for C^, but it cannot be 
presumed that the two functions have other features in 
common . 

Cg and Cg are a measure of the positive and of the 
negative part, respectively, of the source and sink terms in 
the transport equation of the turbulence length scale. The 
former represents the growth of £. as a result of the dis- 
sipation of smaller eddies, the latter represents the 
reduction of Z from the breaking of large eddies . 

C„ and C„ are taken to be functions of R. having 
So tn 

similar behavior as and Cg which take on a constant 
asymptotic value as R^ becomes large. More details of the 
behavior of the coefficient functions will be discussed in 
Chapter IV. 
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6 . 


BOUNDARY CONDITIONS 


It has been mentioned that the governing equations 
(Equations 10, 11, 18, and 21) can be expressed as a general 
elliptical equation having the form of Equation 13. Solu- 
tions of the elliptical equation require specifying the 
boundary condition on a complete contour enclosing the 
region in which a solution is desired. The problem of 
interest is thus a boundary-value problem and the prescribed 
boundary conditions used in this study are discussed below. 

Inlet Boundary Conditions 

The inlet is considered sufficiently far upstream 
that the flow is undisturbed by the building. Thus the 
inlet velocity is assumed logarithmic. Equation 1, in 
keeping with the wind profile of a neutral atmospheric 
boundary layer over flat homogeneous terrain as described in 
Section I. The dependent variable to be specified is the 
stream function which can be determined from 


= J~ pu(z) dz 

o 

Assuming the no-slip condition at z = 0 and carrying out the 
integration gives 


Uj 


^ = P IT 


(z + z,) 


jn - 1 

Z 0 


+ z, 


( 25 ) 


The vorticity, oo, is given by Equation 9, which can be 
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rewritten as: 


(0 = 


p[ax2 9Z2^ 


( 26 ) 


The vorticity at the inlet boundary is thus calculated 
according to Equation 26. 

The boundary condition for the turbulent length 
scale is based on a mixing-length hypothesis which has been 
demonstrated as valid for boundary layer flows [21] such as 
is assumed to exist at the inlet boundary. In the neighbor- 
hood of the wall the mixing length is generally assumed to 
be linearly related to the normal distance from the surface. 
For a very rough surface, the mixing length hypothesis does 
not go to zero but approaches a size on the order of the 
roughness length, zq. The following relationship is assumed 


£ = k (z + z 0 ) 


(27) 


The boundary condition for the turbulence kinetic energy, K, 
is derived from the assumption of constant shear layer. 

T. = T = pu* = constant (28) 

t w ^ 

In terms of Prandtl-Kolmogorov formula. Equation 17, 


= C pK^"'^ £ 
t 


^3u' 


3z 


(29) 


together with Equations 1, 27, 28, and 29 gives 


K = 


C 

'• y-* 


(30) 
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Outlet Boundary Conditions 


Numerous boundary conditions for ijj and to at the out- 
flow boundary condition have been investigated [23] . For 
the flow for downstream from the surface obstacle, it is 
expected that 


3x 


0 


(31) 


Equation 31 implies 


9x2 


0 


(32) 


and hence. 


to 


1 

p 9z2 


(33) 


If = 0 along the lower boundary and if either ip or 
is fixed along the upper boundary, then the setting of 
= 0 contradicts Equation 33 if to 0, or simply restates 
Equations 33 if to = 0 at the upper boundary. 

The gradient of vorticity in the flow direction is 
assumed to vanish at the outlet boundary. 


9(0 

9x 


0 


(34) 


The corresponding conditions for K and IL at the outflow 
boundary are assumed to be 


9K 

9x 


0 


(35) 
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and 


3Jl 

- 


0 


(36) 


which describes an equilibrium boundary layer at the outlet 
boundary. 


Upper Boundary Conditions 

At the free boundary, ip is specified by = w = 0 
and = u = u^, u^ is the inflow velocity at the upper 
boundary, i.e.. 


u 


oo 


k zo 


which assumes the upper boundary is far enough away from the 
obstruction. This boundary condition implies 

ip = constant (37) 

at the upper boundary. 

According to Equation 26, the vorticity on the upper 
boundary thus becomes zero. In the free atmosphere field, 
however, o) = 0 may be unrealistic. A less restrictive con- 
dition for w is therefore assumed: 


3(jo 

3z 


0 


(38) 


If the upper boundary is far from the obstruction so that 
w = 0 on the boundary, = 0 on the boundary and w = 

0 X d Z 

9 ^ u. 

Equation 38 implies = 0. The upper boundary condition 


32 



for ip is 



( 39 ) 


The condition for K is described as a constant on 
the upper boundary which is expected to be neutrally stable. 

K = constant (40) 

Similarly, I is expected to be constant along the upper 
boundary 

S, = constant (41) 

Lower Boundary Conditions 

The solid wall of the region of interest is a stream- 
line along which must be constant. The conventional choice 
is 


ij; = 0 


(42) 


The evaluation of wall vorticity is extremely important and 
can also be problematic, as it is produced by the no-slip 
condition at the wall and thus its specification drives the 
flow. By differentiating the definition of vorticity, one 
obtains 


9z 3x9z 9z2 9x[9zJ p d'z^ 

From the continuity equation. Equation 3, 
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9w 

dz 


3u 

3x 


( 44 ) 


Substituting Equation 44 into Equation 43 gives 

3g) _ 3 _ 3 ^u 

3z ”p 3z^ 3x2 


The last term is zero because of the no-slip conditions, 
hence. 


iiii = ih 

3z ~ p 3z 3 


(46) 


The treatment of the vorticity condition at the upper 
corners of the obstacle requires special attention and will 
be discussed in Chapter III. 

It was assumed that K obeys Equation 30 along the 
wall. By differentiating Equation 1 with respect to Z 

+ ^»>|^ 

Furthermore, at the wall, z = 0 and = - to, thus 

o Z 

u^ = kzo (-to) 

substituting into Equation 30 gives 


^wall 



(47) 


The length scale on the lower boundary was prescribed by 


Jl = kz 0 


(48) 
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CHAPTER III 


NUMERICAL ANALYSIS 

The governing equations together with the boundary 
conditions and transport equations of turbulence kinetic 
energy and turbulence length scale constitute a closed set 
of nonlinear partial differential equations, given by 
Equations 10, 11, 18 and 21. This set of equations has been 
solved by utilizing and extending the numerical procedure of 
Gosman, et al. [11] , which basically develops a set of 
algebraic finite-difference equations by integrating Equa- 
tions 10, 11, 18 and 21 over a finite element of area. 
Because of the nonlinear character of the resulting finite 
difference equations, they are solved by an iterative, 
successive substitution technique. For details of the 
derivation of the finite difference equations, the reader is 
referred to References [11] and [14] . 

1. CONTROL VOLUME AND COORDINATE SYSTEM 

The physical and numerical coordinate system for the 
present study is chosen such that the origin is located at 
the lower left corner of the block. Figure 1, page 3. The 
x-axis is chosen positive in the downstream direction par- 
allel to the surface and the z-axis is perpendicular to the 
x-axis, chosen positive upward and aligned with the left 
face of the block. Figure 1. The subscripts I and J are 
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used to index the x-direction and z-direction respectively 
in the numerical calculation procedure. In the iteration 
process the field was swept from the upper boundary to the 
lower boundary beginning at the inlet boundary and proceeding 
in the increasing I-direction. 

The distribution of the grid points is shown in 
Figure 8. As indicated, a variable mesh was used which 
gradually decreased in size near the wall and in the 
vicinity of the step. 

2, BOUNDARY CONDITIONS TREATED IN THE 
NUMERICAL SOLUTION 


The boundary conditions have already been established 
in Chapter II. The transformations of these into the finite 
difference form used in the numerical procedure are described 
next. 

Inlet Flow ; 

Equations 25, 27 and 30 are directly applicable in 

their original form. The vorticity at the inlet boundary is 

. 9 ^w 

treated by the assumption that = 0. Setting 



1 9 ^ 1 ^ 
p 9x2 


2, J 


(49) 


Together with Equations 1, 9 and 26 


1 9_^ _ u* 1 

p 9x2 2 j ~ ^ z + Zq 


(50) 
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where 


1 

p 3x^ 


2,J 


1 

P 




'P 


1,J 


(Ax) 2 



Outlet Flow ; 

Equation 32 is used for the outlet boundary condi- 
tions far downstream from the block. It is rewritten 


9x2 


IN, J 


9x2 


^IN,J 


^’^IN-1,J '^IN-2,J ^ Q 


IN-1, J 


(Ax) 


where x is uniform at the outlet boundary, which yields 


'^IN,J ^'^IN-1,J '^IN-2,J 


(51) 


The finite difference forms of the vorticity, the turbulence 
kinetic energy and the turbulence length scale are written 
in the following: 


“lN,J ^ ‘*^IN-1,J 
^IN,J ^IN-1,J 
^IN,J ^ ^IN-1,J 


(52) 


Upper Boundary : 

The partial differential form of Equations 38, 39, 
40 and 41 are transferred as 


“l,JN - "'l,JN-l 
'^I,JN ^ ’^1,JN 
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( 53 ) 


^I,JN ^1,JN 

9 =9 

I,JN 1,JN 


Lower Boundary ; 

The vorticity at the lower boundary is developed 
from a Taylor-series expansion of the stream function around 
a near wall point (NP) located a distance n from the wall. 

In terms of the wall point (P) , becomes: 





An + 
P 


1 8^ 

2 8n2 


(An) 2 
P 


+ 


1 

6 3n3 


(An) ^ 
P 


+ H.O.T. 

(54) 


By the no-slip condition/ = 0 and ~ ^nd together 

with Equation 46, Equation 54 becomes 


- . - *P> 

“p (An) 2 . P 



(55) 


The vorticity gradient at the wall was approximated by 


f9“l = *^NP ~ ^P 

^9njp An 


(56) 


which yields upon substitution into Equation 54 

_ _ 3 ^ _ 

“P (An) 2 • . p 2 “np 


(57) 


A special treatment of the vorticity at the upper 
corners of the rectangular block is discussed in [13] . The 
back corner is treated similarly as shown below. Based on 
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the idea that separation occurs tangentially to the upstream 
wall, one can apply Equation 56 to the upstream wall of the 
block and obtain the vorticity at the upper front corner, 
(see Figure 9) as: 

— 1 — 

“c p (Ax ) 2 " 2 “w 

and the vorticity at the upper rear corner, as: 

1 - 

” p (Az) 2 “ "2 “n 

3-. SELECTION OF THE COEFFICIENT FUNCTIONS 

Based on the discussion of Section 5, Chapter II, 

the coefficient functions, C , C^., C„, C„, are all expected 

y jj o o 

to take on constant asymptotic values as becomes very 
large, say C^^, Cj^^, Cg^ and Cg^, respectively. Research 
using the two-equation modeling of turbulence has not 
resolved the exact values of those constants which are 
assigned radically different values by different researchers 
depending on the flow situation being analyzed. 

Wolfshtein [22] , modeling a jet impinging normally 
on a wall, found values of C = 0.22, C^ = 0.416 and 

yoo ' Doo 

CTr gff = 1.53. Rodi and Spalding [24], using production of 
turbulence kinetic energy and length scale, KSL , as a depen- 
dent variable instead of I (K - Kil model) , obtained Cg^ = 

0.09 for free shear flows while was taken as 2.0. Ng 

and Spalding [25] compared boundary layer flows near smooth 
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wall predicted with the (K - K£ model) against experimental 
results for homogeneous shear flows in local equilibrium and 


obtained = 0.1, Cy^ = 1.0 with = 1.0. Launder, 

et al. [20] , employing turbulence kinetic energy dissipation 
rate, e = as the dependent variable instead of il 

found C,| = 0.09 and = 1.0 for free turbulent shear 

flows. In view of the fact that no consistent values of the 
coefficient have been reported for use in the present study 
and that experimental data is scarce, it was necessary to 
develop some understanding of the coefficient functions. 

This is discussed in Chapter IV. Throughout this study 


'K 


and are simply assumed to be unity which is 


consistent with most turbulent flow analysis. 


4. ACCURACY, CONVERGENCE AND ECONOMY 


In carrying out the numerical solution, a balance is 
required between the convergence and accuracy of the results 
and the amount of computing time necessary to meet these 
conditions. The accuracy and economy of the solution pro- 
cedure are opposed to one another, i.e., the finer the grid 
size, the more computing time required but the better the 
accuracy or truncation error. 

A very fine grid size is required, especially near 
the wall, where steep gradients occur. To assure accuracy, 
Wolfshtein [22] found by comparing the exact and the finite- 
difference solutions for the cases of Couette flow, impinging 
flow and uniform velocity flow that a finer grid away from 
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the wall did not improve the accuracy of the vorticity 
solution, but that the stream function was sensitive to the 
grid size. In general, coarser grids were used at a 
reasonable distance away from the walls to cut down the com- 
puting time; however, one of the most important factors 
affecting the convergence is the grid spacing. 

Gosman, et al. [11], suggested that the ratio of the 
consecutive intervals between the nodes normal to the wall 
should be as close to unity as possible, especially near the 
wall but in no case should they be less than 1.5. This 
suggestion was adopted in References [13] and [14] , and is 
also adopted in this study. It was found, in keeping with 
Gosman, et al. [11] , that convergence was impaired when a 
grid having a ratio greater than 2 between the node intervals 
was used, especially in the z-direction [14] . Maintaining a 
uniform grid spacing near the wall assured convergence as 
reported in Reference [13] and as experienced in this study. 

A commonly employed remedy against divergence of the 
iteration process is known as under-relaxation. For the 
present problem there was no need for under-relaxation; 
however, in fact, the stream function, related to the vor- 
ticity by a Poisson type equation (Equation 10) , was over- 
relaxed and the convergence of the solution accelerated. 

5. TERMINATION OF COMPUTATION 

The computation was assumed to have converged to a 
sufficiently exact solution of the difference of the 
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dependent variable (p at all points between successive itera- 


tions became small, i.e.. 


N ,N-1 

p - * 



< 


max 


0.05 


The superscript N denotes the nth iteration. (J>p was chosen 
in the denominator as a scaling factor in order to assure 
relatively good convergence in all areas. 
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CHAPTER IV 


RESULTS AND DISCUSSIONS 

1. INVESTIGATION OF COEFFICIENT CONSTANTS 

The influence of the variation of the constants on 
the flow field is studied in the present investigation. 
Taking the ratio of the production term to the dissipation 
term in Equations 20 and 23, one obtains 


K, 


rC -1 

n 2 

f8u. 9u. 

-1 9u. 1 

poo 

% 

1 + 3 


Cr, 

K 

9x. 9x. 

3 1 

J9x. 

£)ooJ 



'D 


'So- 


le C„ 


K 


1 

f9u . 
X 

9u. 

9x : 
3 

9x. 

X 

"C" 

and "1 


9u. 


(60) 


(61) 


where subscripts "C" and "D" denote the production term and 
dissipation term in Equations 20 and 23, respectively, 

C,. , C„ and C„ are discussed in Section 5 in Chapter II, 
and in Section 3 in Chapter III. Thus 


K, 


'S« 


"d 




(62) 


In view of the fact that the last terms of Equations 60 and 


61 are always positive and the assumption that the influence 

^C ^C 

of those last terms on the and — ratios are always xn 

^D ^D 

the same order, it is expected that these ratios should lead 
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to some understanding of the turbulent levels. For example: 

C 

yoo 

1. Decreasing the ratio decreases the produc- 

'^Doo 

tion of turbulent kinetic energy either through 

less production or faster dissipation. 

^Sco 

2. Decreasing the ratio ^ — = — with C ^ constant 
decreases the growth of the length scale. 

The individual behavior of the coefficient constants can be 
described as listed below: 

a. Decreasing decreases turbulent kinetic 

energy production and also increases the growth 
rate of the turbulent length scale. 

b. Increasing reduces the level of turbulent 

kinetic energy by increasing the rate at which 
it is dissipated. 

c. Decreasing Cg^ or increasing decreases the 

rate of growth of the turbulent length scale. 

The influence of C on the turbulent level is 

yoo 

illustrated from Equation 30. The inlet condition of turbu- 
lence kinetic energy is strongly dependent on and thus 
the variation of C will influence the turbulence kinetic 

yoo 

energy over the entire flow field. This is reasonable in 
that the turbulence level in the incident flow should 
influence the flow in the region of interest. The smaller 
C value will cause larger turbulence kinetic energy in the 

yoo 

flow field. This fact appears opposite to the prediction 
"a," above; however, it suggests that the higher the incident 
turbulence the lower the production rate which is not 
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implausible. The behavior of still requires further 

investigation. 

Frost and Harper [26] point out that for the 
atmospheric boundary layer it has been experimentally deter- 
mined that 

a = 2.5 u* 
u ^ 

= 1.9 u* 

for a neutral atmosphere, where 


r zs 

(U- 

“) 

1/2 

U 


r zr 

(^ 

^) 

1/2 

V 


w ~ 

(w* 

“) 

1/2 


These relationships together with Equation 15 give 

K = 5.78 u^ (63) 

Considering the inlet boundary condition for turbulent 
kinetic energy given by Equation 30, one obtains 

C = 0.416 (64) 

yoo 

for the neutrally stable atmosphere. This value of is 

adopted throughout the remainder of this investigation. Two 
cases were chosen to investigate the influence of the 
coefficient constants. They are: 
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I 


Case I : 


C 


yoo 


0.416 


'S«> 


1.0 


C 


Doo 


1.0 




1.0 


Case II; 


C = 0.416 

yoo 


^D«> 0*43 


Cg. = 0-60 




1.40 


yoo gcx3 

The behavior of the ratios, and ^ ^ — , are 

considered rather than the individual variation of the 
coefficients. Figure 10 indicates the turbulent kinetic 
energy of Case II is larger than that of Case I at the 
corresponding position in the flow field. This is expected 
because 


c 


fc 1 

yoo 

> 

yoo 

^Doo 



II 



meaning the ratio of production to dissipation is greater 
for Case II than Case I. Figure 11 shows the turbulent 
scale for Case I is larger than that of Case II at the 
corresponding point in the flow field. This effect is 
caused by the fact 


1 

> 

✓ \ 

yoo goo^ 

I 

‘“yoo^Boo 


It is observed that the qualitative behavior of the 
results shown in Figures 10 and 11 are very similar despite 


48 




49 
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the choice of the coefficient constants but that the 
quantitative results depend on the selection of these 
coefficient constants. 

To establish some physical explanation of the ratio 

Cg^ 

— p; — one can argue that the growth of the turbulent length 
scale is a result of dissipation of energy which creates the 
biggest fatality rate for small eddies. Thus when dissipa- 
tion exceeds production the length scale, which is taken to 
represent some effective mean of the distribution of length 
scales, will increase with the depletion of the smaller 
eddies. The reduction in the length scale, in turn, may be 
thought of as being caused by the tendency of shear stress 
to rupture the large eddies. Therefore, when the production 
of turbulence kinetic energy exceeds dissipations the length 
scale will decrease and when the dissipation exceeds pro- 
duction the length scale will grow. Hence one might suspect 



( 65 ) 


and we assume here the relationships 




1.0 


( 66 ) 


This assumption is at least valid for quantitative study. 
The relationships of coefficient constants are generally 
derived under the consideration of an equilibrium boundary 
layer. Appendix A, but this approach is not suitable in the 
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more complex flow problem presented herein. Due to the 
experiences with the present computer code reported in 
References [13, 14], the program is known to work well when 
all the coefficient constants are chosen to be unity. 

Taking cognizance of these experiences and considering 
Equation 66, the following values are chosen for the 
qualitative studies. 

Case III; 


0.416 


0.416 



The results of a qualitative study with these values and for 
u* = 36.63 cm/sec and = 0.7 cm are reported in the next 
section. 


2. VELOCITY DISTRIBUTION 

The velocity distribution of the flow around a bluff 
obstruction is one of the most interesting topics in current 
numerical fluid mechanics research. Figure 12 shows the 
computed velocity profiles at selected X stations. In the 
region close to the wall, the flow is decelerated as it 
approaches the obstruction and is accelerated as it passes 
over the obstruction. A recirculating flow region which 
extends approximately 11 step heights downstream is created 
behind the obstruction. In the region above the recircula- 
tion zone, the flow is accelerated because of the displace- 
ment of the flow. The flow reattaches near X = 12.25 H and 
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H 

(H = 3.2 m) 

Figure 12. Velocity profiles around an obstruction on the surface. 


the logarithmic boundary layer begins to reestablish down- 
stream. Figure 13, which is a plot of the local velocity 
profiles nondimensionalized with the upstream value clearly 
demonstrates the phenomena described above. The results 
in Figure 13 clearly support the often reported observation 
[12] that the flow tends to equilibrium at a much slower 
rate than the rate at which the disturbance process takes 
place. Figure 12 indicates that at 20 step heights behind 
the obstruction the flow has not redeveloped to the initial 
state. This is supported by Reference [5] which reported 
that 44 step heights were required for the flow to return 
to equilibrium. 


3. TURBULENCE KINETIC ENERGY 

Figure 14 gives some computed information about the 
tendency of turbulent kinetic energy in the flow field. In 
this study, ijj , w, K and H are chosen as the four dependent 
variables. The behavior of any one of these dependent 
variables influences all of the others. Comparing Figures 
14 and 13, one concludes that the deceleration of flow in 
the flow direction tends to increase the turbulent kinetic 
energy and that the acceleration of the flow in the downwind 
direction reduces the local turbulent kinetic energy. The 
extremely small value of the turbulent kinetic energy found 
below z = H at 0.25 H behind the obstruction is believed to 
be caused by the rear face of the obstruction restricting 
the fluctuation in that region. Turbulent kinetic energy 
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Figure 14. The ratio of 
inlet condition at t 




near a solid surface has typically low values because the 
random motion of fluid particles is restricted. 

One can trace the shear layer separating from the 
front face of the block and diffusing outward as it moves 
downstream. This shear layer is highly turbulent as the 
computation clearly portrays. One can also well observe the 
diffusive action of the turbulence kinetic energy in the 
shear layer originating from the corner of the block. 

Figure 14 shows that the turbulent kinetic energy reaches 
its peak value over the reattachment point and then begins 
to decrease and to tend to a stable state downstream. Tur- 
bulent kinetic energy decrease in the flow direction behind 
the reattachment point is almost the same as that of turbu- 
lent kinetic energy increase in front of the reattachment 
point. The turbulence intensity is still considerably in 
excess of the initial value at 17 step heights downstream. 

4. TURBULENCE LENGTH SCALE 

The turbulent length scale employed herein is the 
integral of the two point velocity correlation in the flow 
direction. This is called the integral length scale and is 
considered a measure of the largest eddy size. The turbu- 
lent scale growth was suggested as representing the effect 
of preferential dissipation of ^mailer eddies by turbulence 
kinetic energy dissipation. Considering the ratio of 
turbulence length scale to the initial value. Figure 15, and 
that of velocity to the inlet value. Figure 13, page 55, one 
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can observe that the acceleration in the flow direction 
always shrinks the eddy length. The effect of deceleration 
gives the opposite effect. 

The large length scale near the lower corners of the 
obstruction may be explained in terms of the two vortices 
near the front and rear lower corners, respectively. 

Further discussion of these vortices are given in the next 
section. 


5. VORTICITY AND STREAMLINE PATTERNS 

Both computed vorticity contours. Figure 16, and 
streamline patterns, Figure 17, indicate that a small down- 
wash zone occurs near the front lower corner and a large 
recirculation zone occurs behind the obstruction. For 
theoretical analysis the boundaries of both zones are chosen 
as the streamline, = 0. The zero streamlines divide the 
mixing region, which comprises the original and part of the 
new shear layer, and the outer flow from the downwash zone 
and recirculation zone where the fluid recirculates as a 
large eddy . 

6. COMPARISON OF COMPUTED RESULTS 
WITH EXPERIMENTAL DATA 

Two data sets obtained from the field tests (see 
Chapter I) , data sets #8505 and #8512 [5] are chosen for 
comparison with the calculated results. Figure 18 shows non- 
dimensional velocity profiles, the local velocity divided by 
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: Computed Results 

O ; Data set #8505 (u^^ measured at tower #1 [5]) 

□ : Data set #8512 (uj^^ measured at tower #1 [5]) 

• : Data set #8505 (uj^^ measured at tower #6 [5]) 

■ : Data set #8512 (Uj^^ is measured at tower #6 [5]) 

Figure 18. Comparison of the ratio u/uj^j^» 



the inlet velocity at the same altitude, and Figure 19 
shows nondimensional turbulence kinetic energy profiles, 
the local value divided by the inlet magnitude at the same 
altitude at several x stations down the array. Observing 
Figure 18, one notes that the computed results show reason- 
able agreement at the upstream stations but depart 
appreciably at the downstream locations. Figure 18 indi- 
cates the experimental velocity deficit, u^^ - u, decays 
by 5 block heights downwind from the block. This value 
is far smaller than that obtained in the wind tunnel studies 
(see Chapter I) and of that found in Reference [5] which 
conducted smoke tests and showed that reattachment occurs 
at 12 H. The computed results indicate that the velocity 
defect decays in the same manner as the observation of the 
wind tunnel tests [28] . The details of this comparison will 
be performed later. 

At the location x = 17.75 H, the nondimensional 

velocity has an unexpected shift from the expected value of 

unity at the lower levels. In the field tests, there is a 

small hill in front of Tower #5 which is located at 16.5 H 

behind the block. This hill is expected to create an 

acceleration due to mass conservation and the velocity u 

will thus exceed the upstream value causing ■ to be greater 

^in 

than unity. If one assumes the measurement at Tower #6, 
Figure 4, page 13, as the reference, one can observe from 
Figure 19 that the computed results are always larger than 
the results from the field tests and the acceleration 
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in 

: Computed Results 

O : Data set #8505 [5] 

□ : Data set #8512 [5] 


Figure 19. Comparison of the ratio K/K.. 


phenomena disappears. The explanation is that the presence 

of the trees ahead of the block at the field site create a 

velocity deficit. The velocity measured at Tower #6 [5] is 

therefore larger than that measured at Tower #1 [5] due to 

this wake from the trees, and thus the ratio where u. 

'^in 

is measured at Tower #6, will always be lower than the ratio 

where u. is determined from Tower #1. 
xn 

Figure 19 shows that the computed turbulence kinetic 
energy agrees well with the experimental data at locations 
above four block heights, but does not demonstrate the 
regions of high turbulence found from computation and v/ind 
tunnel experiments. The experimental data (Figure 19) 
implies that the flow field is always in a local equilibrium 
state even in the presence of the obstruction which is 
contradictory to the observation of the wind tunnel tests 
[28, 29] and to the current results. This may result from 
the fact that the lateral scale of turbulence in the 
atmosphere is much greater than that simulated in the wind 
tunnel. In the numerical model lateral turbulence is, of 
course, neglected. The rms value of the wind direction 
measured in the field studies [5] is on the order of ±56 
degrees. This is expected to cause some disagreements 
between the field test data and the two-dimensional computed 
results. 

Figure 20 provides a comparison of the value u'w' of 
the current computed results (calculated by u'w' = 

\>u.(-|^ + 1^)) with the results of the field study. In 

^ o Z d X 
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: Computed Results 

O : Data set #8505 [5] 

a : Data set #8512 [5] 


Figure 20. Comparison of the ratio (u'w')/(u*). 
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Figure 20 the values of u'w' have been nondimensionalized by 
dividing by the square of friction velocity, — The 
agreement is very poor. The presence of the trees, which 
are approximately 3.5 H high and approximately 40 H upstream 
of the block in the actual field site studied, create a 
mixing region and are thus believed to produce high values 
of u'w' at heights near that of the tree tops. The current 
computed results show that the value of v\j. tends to 

a constant far downstream from the block. This coincides 
with the concept of a constant shear stress layer for 
neutrally stable atmospheres (see Chapter II, Section 2 ). 

As mentioned earlier, wind tunnel tests have been 
conducted which also provide data for comparison with the 
current computed results. Figure 21 indicates the decay of 
the mean velocity defect along the wake centerline as 
computed and as measured in the wind tunnel [28] . Current 
results coincide almost perfectly with the experimental 
data. Figure 21 shows that two-dimensional wakes, associ- 
ated with a long building, having the ratio height/width/ 
depth = 1/8/1, is larger than the three-dimensional wakes 
associated with short building, H/W/D = 1/2/1. The very 
good agreement of the computed results with the experimental 
data for the long block geometry tends to confirm two- 
dimensionality along the centerline for large aspect, ratios. 

A comparison of the vertical profiles of the mean 
velocity defect at selected ^ stations are shown in Figures 
22 and 23. These indicate that both the computed results 
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■ : Cermak, et al. [28] (nearly two-dimensional model) 

: Current results 

: Cermak, et al. [28] (three-dimensional model) 


Figure 21. Decay of mean velocity defect along the wake center line 
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Figure 23. Vertical profiles 
two-dimensional model) [28 



and the experimental data display a peak value at approxi- 
mately the height 0.8 H in the vicinity of the block. 

However, the shear layer appears to diffuse outward slower 
in the computed results than that in wind tunnel studies. 

Also the experimental data appear to recover slower than the 
prediction of this study. The reason is that the flow field, 
which is considered in the computer study, is specified at 
20 H downstream from the block to have zero vertical 
velocity. Thus, diffusion may be curtailed by the boundary 
condition. 

The magnitude of the measured and computed velocity 
defects are approximately a factor of two different. How- 
ever, the data of Figure 21 taken from the same reference 
for apparently the same test condition do not have this 
factor of two disagreement. It is therefore believed that 
the scale given in [28] for Figure 23 is off by a factor of 
two. If this is the case the computed results show excellent 
agreement with the experimental results. 

Other factors which are expected to cause some dis- 
agreement between the wind tunnel data and the theoretical 
results but not as large as a factor of two difference in 
magnitude are as follows. The logarithmic velocity profile, 
which was used as an inlet condition in this study, does not 
coincide exactly with the power profile used in the wind 
tunnel study. A comparison of these profiles is given in 
Figure 24, where the inlet velocity at six block heights is 
chosen as the reference velocity. The inlet turbulence 
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Figure 24. Comparison of logarithmic velocity profile with 
the power law velocity profile. 
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kinetic energy is based on the concept of a constant shear 
layer for a neutrally stable atmospheric boundary layer and 
thus the gradients of velocity affect the value of the inlet 
turbulence kinetic energy (Equations 28 and 29) . Before 
discussing this further, the influence of different inlet 
turbulence kinetic energy levels are discussed. For this 
purpose, two cases are chosen — Case III and Case IV which 
assxames C = C„ =1.0 and C_ = C„ =1.0. The effects of 

JJOO IJOO goo goo 

increased inlet turbulence kinetic energy can be assessed 
from Figure 25 which compares the velocity computed for 
Case III with that for Case IV. From Equation 30, the inlet 
turbulence kinetic energy for Case III, is larger than 

that for Case IV, Inspection of Figure 25 shows that 

the separation bubble for Case III is slightly larger than 
that for Case IV. One can conclude that the larger the 
value of the inlet turbulence kinetic energy the larger the 
separation bubble. This is not consistent with intuition 
since large turbulence level would be expected to force 

AK^ 

reattachment sooner. However, it seems that the ratio ^ — , 

^in 

where is the local turbulence kinetic energy and is 

the inlet value (see Appendix B) , must be considered. 

Observing Figure 25, one further notes that the velocity 
Uin ” u 

defect, — of Case III, is smaller than that of Case IV 

in the upper region and larger, in the lower region. This 
trend is also believed to be associated with the variation 
in inlet turbulence kinetic energy. In Reference [28] the 

/^l 2 \ 1/2 

upwind turbulence intensity, r is approximately 8% 
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Figure 25. Comparison of velocity profile in the wake region 




which is smaller than the theoretical case investigated in 
this study, - — = 12.5%. 

Uoo 

If the turbulence were isotropic, = (u'^)^'^^. 

However, it is expected that the fluctuation velocity in the 
flow direction is larger than that in the vertical direction 
so that u ' 2 > w' 2 and thus 


u„ 


(K) 


1/2 


u 


For example, if w^ = 0.5 u'2, then 


(K) 


1/2 ^ 



+ 





1/2 


and 




1.5 


1/2 


(K) 


1/2 


(K)1/2 


It is therefore to be expected that the velocity defect, 

Ufn - u 

, of the wind tunnel studies [28] is larger in the 

Uoo 

upper region and smaller in the lower region from those pre- 
dicted by the numerical computation. 

Hunt [8] develops a theory for a two-dimensional 
wake behind a surface block. A comparison of the current 
results with Hunt's theory is shown in Figure 26. Reference 
[8] compares the theory with the experimental data from 
Counihan [8] which is obtained for a power profile with 
m = 0.125 and shows the data to collapse onto one similarity 
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Hunt's theory [8] 

Current results (x/H = 20) 
Current results (x/H = 7.0) 
Current results along the 
center line of mixing layer 


Figure 26. Comparison of Hunt's theory with current 
computed results (m = 0.125). 





curve. The computed results here, however, do not show the 
same tendency to collapse onto a single curve as suggested 
by Hunt's theory. Reference [28] also found that their data 
did not fall on a similarity curve. It should be noted that 
the computed results along the approximate center line of 
the mixing layer at z* 1.25 H agree relatively well with 
Hunt's theory. 

Reference [ 8 ], whose data agree with Hunt's theory 
well , shows a separation bubble which extends to x/H = 6.0 
and is considerably smaller than that of the computed 
results for which the reattachment point occurs at x/H = 
11.5. This difference is believed to be a reason why the 
computed results reported in this study depart from Hunt's 
theory. 

It is useful to investigate the character of the 
computed shear layer shed from the leading edge of the 
building. Many experiments [14, 29, 30] have shown the 
velocity in the shear layer to be governed by the erf 
function profile. In this study, the upper boundary of the 
mixing region is defined as u/u^ = 1 . 0 , where u^ is taken as 

the inlet velocity at the upper boundary of the region 

U 6 H + Zri u 1 

(i.e., u„ = in -) . The line — = - 5 - as conven- 

k Zq *^00 ^ 

tionally done, treated as the jet axis of the mixing layer 
theory, — = -^( 1 +erf ?) where 5 = o y is measured from 

jet- axis and a is a constant to be determined. Figure 27 
indicates the jet-axis and upper boundary of mixing layer. 
The wedge fence model [29] gives a higher upper mixing 
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— : Current results 

— ; Upper jet boundary of wedge fence model [29] 
— : Jet axis of wedge fence model [29] 


Figure 27. Upper mixing boundary and jet-axis. 




boundary and a higher jet- axis than the computed results for 
the rectangular block. This is reasonable in that the wake 
region of the fence extends much higher than that of a 
forward facing step and the rectangular block represents an 
intermediate case. The height of the wake region is 
obviously a function of the recirculating flow and resulting 
pressure build-up for the different geometries. It should 
be noted that the current computed results of a wake 
approximately 2 H in height coincides with the results of 
the smoke trace study carried out in Reference [5] . 

A comparison of the mixing layer theory with the 
current computed results and the values of the computed 
spreading parameter, a, are shown in Figure 28. Figure 29 
is a plot of the spreading parameter a with respect to x/H. 
The results from the wedge fence model [29] is also shown in 
Figure 29. The explanation of the difference between the 
computed results and the experimental data could be that the 
value of a is dependent on the geometry of the bluff body 
from which separation occurs. Evidence that a depends on 
geometry is given in Reference [30] . Another possible 
reason for the difference between the computed a and the 
measured values of [29] is that the over prediction of the 
turbulence kinetic energy in the flow field may force the 
flow to bend to the ground early and make the mixing region 
smaller than that observed in Reference [29] . 

Figure 30 shows the computed turbulent properties 

— ~ and (or v.. (— + — ) ) . It is observed that the peak 

^ 9z ax 
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Figure 28. Comparison of mixing layer theory and current computed 
results. 
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Figure 29. Variation of similarity parameter with the 
ratio x/H. 




value turbulence kinetic energy occurs along the axis of 
mixing layer and then approaches a constant outside the 
upper mixing boundary. The former phenomena coincides with 
the observation in Reference [29] for the turbulence 
intensity which is proportional to the turbulence kinetic 
energy (see Figure 31 ) , but the turbulence intensity of [29] 
approaches zero rather than a constant as in Figure 30. 

The explanation is that the flow field considered in this 
study is a neutrally stable atmospheric layer which possesses 
constant shear stress in the undisturbed case and thus has 
constant turbulence kinetic energy. 

The turbulent shear stress in Figure 30 has a 
different sign from that in Figure 31. This is simply 
because of the sign convention chosen for displaying the 
data. 


7. CONCLUSIONS 

The present approach of solving the ensemble 
averaged two-dimensional incompressible turbulent Navier- 
Stokes equations as applicable to a neutrally stable 
atmospheric flow over a two-dimensional block has been shown 
to yield results in agreement with the experimental findings. 
The two equation model (K - £ model) of turbulence has been 
demonstrated to produce realistic results in the recircula- 
tion flow field reported in this study. 

The influence of the coefficients' ratio on the 
turbulence kinetic energy and turbulence length scale of 
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the flow field of interest and on the predicted results have 
been computationally investigated, and a set of constants 


has been determined for the problem under study. 

C 

qoo 

results of the study show that the ratio ^ — 


close to unity. 


The 

has a value 
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APPENDIX A 


COEFFICIENT CONSTANTS OBTAINED FROM 
BOUNDARY LAYER APPROACH 


Considering the transport equations of turbulence 
kinetic energy and turbulence length scale, e.g.. Equations 
18 and 21, and taking the boundary layer approach. Equations 
18 and 21 are reduced to 
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Introducing the concept of a constant shear stress 
layer such that turbulence kinetic energy is constant and of 
the turbulence Reynolds number, R^ = ^ Tjj, is very large, 

'^yoo ^ 

Equation 24. Equation A-1 simplifies to 


't[3zj 


= ^dP 
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(A-3) 


Thus the production and dissipation of turbulence kinetic 
energy are in equilibrium. 

Wind profile in surface layer is described as 
Equation 1, such that 
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where Jl = k(z + z^) = 0.4(z + z^) . 

Substitute into Equation A- 3 
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together with Equation 30 
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In the case of negligible convection effect, e.g., vicinity 
to the wall. Equation A-2 is reduced to 
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so that 

+ Cg - Cp • Cg = 0 (A-8) 


or 
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(A- 9) 


On the other hand, one can find a relationship 
between coefficient constants in equilibrium state from 
Equation A-7. That is 


CgK^/^p ^ 
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or 
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(A-11) 



Based on the consideration of boundary layer approach, one 
has already obtained: 

1. Equation A-6 together with Equation A-9 is valid 
for the case in which convection effect can be 
neglected. 

2. Equation A-6 together with Equation A-11 is 
available for the equilibrium state. 

It is necessary to mention that those results are 

3 u 

obtained by the assumption — = 0 so that it is not accept- 

o X 

able in a complicated flow situation, e.g., recirculation 
flow. The turbulence length scale is assumed to obey 
Prandtl mixing length hypothesis vicinity to the wall. 
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APPENDIX B 


PREDICTION OF THE SIZE OF SEPARATION 
BUBBLE BEHIND A BLOCK 

The local turbulence kinetic energy, may be 

denoted by the following: 

"il ,2 = ''1,1 "c - ''d 

where the subscripts "1," ”2,” "C" and "D" denote the up- 
stream condition, the downstream condition, the production 
and the dissipation of the local turbulence kinetic energy, 
respectively. According to Equation B-1, the following 
equation is obtained: 

^g-,2 ~ ^)l,l 


One can define a ratio R by the following: 
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(B-2) 


One can also denote the relationship of the ratio R of two 
different cases, say case I and case II, as follows: 


II 


< Vi _ 
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(B-3) 
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According to dimensional analysis (refer to Equation 20) , 
one denotes 






(B-4) 


Substituting Equation B-4 into Equation B-3, one obtains: 
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(Cp>i(K 
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(B-5) 


Assuming z for the same inlet boundary condition of 


K. 


the turbulence length scale and assiiming the ratio ^ is the 
same for the same geometry array of the flow field. Equation 
B-5 is simplified in the following: 


R, 
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(B-6) 


Based on Equation B-6, the following discussions are given. 
Discussion I : 

If the inlet boundary conditions of two cases are 
the same, it is reasonable to assume the local turbulences 
of these two cases are always in the same order. Basing on 
this assumption. Equation B-6 is simplified as follows. 




(B-7) 
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one expects Rj is larger than 


If Cpj is larger than 
Rjj. Due to the definition of the ratio R, one may argue 
that the larger the ratio R the larger the turbulence 
tendency in the flow field. Basing on this argioment and the 
fact that the laminar flow always creates a larger separa- 
tion bubble than the turbulent one, one concludes that the 
larger the value the smaller the separation bubble for 
the cases which possess the same inlet turbulence kinetic 
energy . 

In this study, it has been discovered that the 
separation bubble extends approximately to 6 H behind the 
block for the case = 1.0 and approximately 11.5 H for the 
case = 0.416. 

Discussion II ; 

Due to the reason that the problem concerned in this 

study is a boundary- value problem (see Chapter II, Section 

6 ) , one may argue that the larger inlet turbulence kinetic 

energy always creates the larger local value at a specified 

location in the flow field. According to Equation 30, one 

expects that the larger value the larger inlet turbulence 

kinetic energy. For the cases which possess the same ratio 
C 

but have different inlet turbulence kinetic energy, one 
^D 

expects that the larger inlet turbulence kinetic energy the 

larger value of Considering Equation B-6, assuming 

K. > K. ,, and basing on the arguments stated above, 

in, II in, I 

one expects that Cj^ ^S,,II ^ ^S,,I thus the 
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magnitude of R to compare with that of is ambiguous 

± c 

(Equation B-6) . In this study the ratio is chosen to be 

unity. The separation bubble extends to 11.5 H behind the 

block for the case Cj^ = 0.416 and extends to approximately 

11.0 H for the case = 1.0. One may conclude that the 

size of the separation bubble does not change very much even 

C 

though the variation of C_ value is large if the ratio is 
kept constant. 
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